In [1]:
import numpy as np
import scipy.io as sio
import matplotlib.pyplot as plt
from numpy.fft import  fft
%matplotlib inline

In [2]:
data = sio.loadmat('sunspotData.mat')
years, sunspots = data['years'], data['sunspots']
start, end = min(years)[0],  max(years)[0]
ave_sunspot = np.empty(end - start + 1)
for i, year in enumerate(range(start, end + 1)):
    ave_sunspot[i] = np.mean(sunspots[years == year])

plt.plot(np.arange(start, end + 1), ave_sunspot)
plt.xlabel('t [year]')
plt.ylabel('sunspots')


Out[2]:
<matplotlib.text.Text at 0x7f51485a5cf8>

We can clearly see in that the number of sunspots is periodic with a period of more or less 10 years


In [3]:
plt.plot(ave_sunspot[:-1], ave_sunspot[1:])
plt.xlabel(r'sunspots$_{t - 1}$')
plt.ylabel(r'sunspots$_{t}$')


Out[3]:
<matplotlib.text.Text at 0x7f51485a5860>

In [18]:
omega= fft(sunspotsdef estimated_autocorrelation(x):
    n = len(x)
    variance = x.var()
    x = x-x.mean()
    r = N.correlate(x, x, mode = 'full')[-n:]
    #assert N.allclose(r, N.array([(x[:n-k]*x[-(n-k):]).sum() for k in range(n)]))
    result = r/(variance*(N.arange(n, 0, -1)))
    return result).real

In [19]:
plt.plot(np.abs(np.real_if_close(omega)))


Out[19]:
[<matplotlib.lines.Line2D at 0x7f51481bdf60>]

In [20]:
def estimated_autocorrelation(x):
    n = len(x)
    variance = x.var()
    r = np.correlate(x, x, mode = 'full')[-n:]
    result = r/(variance*(np.arange(n, 0, -1)))
    return result

In [27]:
acorr = estimated_autocorrelation(ave_sunspot.flatten())
plt.plot(range(len(acorr)), acorr)


Out[27]:
[<matplotlib.lines.Line2D at 0x7f5147e45b38>]

In [ ]: